subroutine usdefnc2( un, u,  ub, uf, ucb, ucsb, ucf, ucsf &
                         , pn,p, pb, pf, pcb, pcsb, pcf, pcsf   &
                         , ron,rho,rhob,rhof,rhocb,rhocsb,rhocf,rhocsf   &
                         ,ds1,ds2,al,ar,dt,dtmin,epsilon2)
implicit none
	!IMPLICIT REAL *8(a-h,o-z)
	integer k,rank
	real *8 al(4,3),ar(4,3)
	real *8 gi(3,3)
real *8	teta1Alfa1l,l1l
real *8	teta2Alfa1l
real *8	teta3Alfa1l
real *8	teta1Alfa1r,l1r
real *8	teta2Alfa1r
real *8	teta3Alfa1r
real *8 char3l,char3r

    real *8 esm(3)
	real *8  un(3), u(3),  ub(3), uf(3), ucb(3), ucsb(3), ucf(3), ucsf(3)
	real *8 pn,p, pb, pf, pcb, pcsb, pcf, pcsf
	real *8 rhon,rho,rhob,rhof,rhocb,rhocsb,rhocf,rhocsf
	real *8 ds1,ds2,dt,dtmin
	real *8 uq(3,10)
    real *8 ei,z1,vtmp,dttmp
    real *8 q,r,qn,rn,u1x,ron,vf,wf,vb,wb
    real *8 char3f,soundf,gf,sf,qf,qcf,qcsf,char2f,qwellf,qmaxf,qminf
    real *8 char3cf
    real *8 char3csf
    real *8 char3b,soundb,gb,sb,rb,rcb,rcsb,char2b,rmaxb,rminb
    real *8 char3cb
    real *8 char3csb
! SK EXTRA VARIABLES
    real *8 qb,qcb,qcsb,rf,rcf,rcsf,rwellb
	real *8 qwellb,qmaxb,qminb,rwellf,rmaxf,rminf
	real *8 gl,gr,rnr,rnl,qnl,qnr,char1f,char1b,snd,sndf,sndb
	real *8 smachl,smachr,smach,gn,rc1,rc2,rc3,rc4,qc1,qc2,qc3,qc4,var1,var2
! END SK
    real *8 charc1,charc2,charc3
    real *8 a(3,3)
    real *8 epsilon
    real *8 es(3)
    real *8 deff,defb
    real *8 dtmint,gam,gam2,gamx,gamy,s,scb,scf,scsb,scsf,sn,sound
    integer i,j,iback,l
    external deff
    external defb
    real *8 epsilon1,epsilon2

    real(8) :: alfa, sn1,sn2, ut11,ut21,ut12,ut13,ut22,ut23, pn1,pn2, ron1,ron2

! SK
!    epsilon=epsilon2
	iback=1
    !if (pcsb.lt.0.95d5) epsilon=0.2
    !if (pcsf.lt.0.95d5) epsilon=0.2
    !if (pcsb.gt.1.15d5) epsilon=0.2
    !if (pcsf.gt.1.15d5) epsilon=0.2
    gam=1.4d0	! Heat capacity ratio for a diatomic gas
    gamx=2.0d0*dsqrt(gam)/(gam-1.0d0)
    gamy=((gam-1.0d0)/(2.0d0*gam))
    gam2=2.0d0*gam
    !es=esm/dsqrt(esm(1)**2+esm(2)**2+esm(3)**2)
!	es(1:3)=al(1:3,1)
!	call coor(es(1),es(2),es(3),a)
! 1-u 2-ub 3-uf 4-ucb 5-ucsb 6-ucf 7-ucsf
!    call  transform(u,ut(1,1),a)
!    call  transform(ub,ut(1,2),a)
!    call  transform(uf,ut(1,3),a)
!    call  transform(ucb,ut(1,4),a)
!    call  transform(ucsb,ut(1,5),a)
!    call  transform(ucf,ut(1,6),a)
!    call  transform(ucsf,ut(1,7),a)

	call invg(al,gi)


	teta1Alfa1l=0d0
	teta2Alfa1l=0d0
	teta3Alfa1l=0d0
	do j=1,3
	teta1Alfa1l=teta1Alfa1l+al(j,1)*al(j,1)
	teta2Alfa1l=teta2Alfa1l+al(j,1)*al(j,2)
	teta3Alfa1l=teta3Alfa1l+al(j,1)*al(j,3)

!	teta1Alfa1l=teta1Alfa1l+(gi(j,1)*al(1,j)+al(j,1)*gi(j,1))/2d0
!	teta2Alfa1l=teta2Alfa1l+(gi(j,1)*al(2,j)+al(j,1)*gi(j,2))/2d0
!	teta3Alfa1l=teta3Alfa1l+(gi(j,1)*al(3,j)+al(j,1)*gi(j,3))/2d0

	end do

	l1l=dsqrt(teta1Alfa1l)
	uq=0d0
! 1-u 2-ub  4-ucb 5-ucsb
	do i=1,3
	do j=1,3
	uq(i,1)=uq(i,1)+al(j,i)*u(j)
	uq(i,2)=uq(i,2)+al(j,i)*ub(j)
	uq(i,3)=uq(i,3)+al(j,i)*ucb(j)
	uq(i,4)=uq(i,4)+al(j,i)*ucsb(j)
	end do
	end do
	do j=1,4
	uq(2,j)=uq(2,j)-teta2Alfa1l/teta1Alfa1l*uq(1,j)
	uq(3,j)=uq(3,j)-teta3Alfa1l/teta1Alfa1l*uq(1,j)
	end do
	call invg(ar,gi)
	teta1Alfa1r=0d0
	teta2Alfa1r=0d0
	teta3Alfa1r=0d0
	do j=1,3
	teta1Alfa1r=teta1Alfa1r+ar(j,1)*ar(j,1)
	teta2Alfa1r=teta2Alfa1r+ar(j,1)*ar(j,2)
	teta3Alfa1r=teta3Alfa1r+ar(j,1)*ar(j,3)
!!!	teta1Alfa1r=teta1Alfa1r+ar(j,1)*gi(1,j)
!!!	teta2Alfa1r=teta2Alfa1r+ar(j,1)*gi(2,j)
!!!	teta3Alfa1r=teta3Alfa1r+ar(j,1)*gi(3,j)

!	teta1Alfa1r=teta1Alfa1r+(gi(j,1)*ar(1,j)+ar(j,1)*gi(1,j))/2d0
!	teta2Alfa1r=teta2Alfa1r+(gi(j,1)*ar(2,j)+ar(j,1)*gi(2,j))/2d0
!	teta3Alfa1r=teta3Alfa1r+(gi(j,1)*ar(3,j)+ar(j,1)*gi(3,j))/2d0

	end do




	l1r=dsqrt(teta1Alfa1r)
	
! 6-u 7-uf  9-ucf 10-ucsf
	do i=1,3
	do j=1,3
	uq(i,5)=uq(i,5)+ar(j,i)*u(j)
	uq(i,6)=uq(i,6)+ar(j,i)*uf(j)
	uq(i,7)=uq(i,7)+ar(j,i)*ucf(j)
	uq(i,8)=uq(i,8)+ar(j,i)*ucsf(j)
	end do
	end do
	do j=5,8
	uq(2,j)=uq(2,j)-teta2Alfa1r/teta1Alfa1r*uq(1,j)
	uq(3,j)=uq(3,j)-teta3Alfa1r/teta1Alfa1r*uq(1,j)
	end do


    char3l =uq(1,1)
	char3r =uq(1,5)
    char3b=uq(1,2)
    char3f=uq(1,6)
    char3cb=uq(1,3) 
    char3csb=uq(1,4)
    char3cf=uq(1,7) 
    char3csf=uq(1,8)
    s=dlog(p)-gam*dlog(rho)
    sb=dlog(pb)-gam*dlog(rhob)
    sf=dlog(pf)-gam*dlog(rhof)
    scb=dlog(pcb)-gam*dlog(rhocb)
    scsb=dlog(pcsb)-gam*dlog(rhocsb)
    scf=dlog(pcf)-gam*dlog(rhocf)
    scsf=dlog(pcsf)-gam*dlog(rhocsf)
	gf=gamx*dexp(scsf/gam2)*l1r
	gb=gamx*dexp(scsb/gam2)*l1l
    soundf=dsqrt(gam*pcsf/rhocsf)*dsqrt(teta1Alfa1r)
    soundb=dsqrt(gam*pcsb/rhocsb)*dsqrt(teta1Alfa1l)

    sndb=dsqrt(gam*pb/rhob)*dsqrt(teta1Alfa1l)
    snd=dsqrt(gam*p/rho)*(dsqrt(teta1Alfa1r)+dsqrt(teta1Alfa1l))/2d0
    sndf=dsqrt(gam*pf/rhof)*dsqrt(teta1Alfa1r)


! SK USE PANIKOVSKI FOR THE CENTRED EXPANSION FAN 
! RIGHT CELL
	epsilon=epsilon2
	if(char3r.le.0 .and. char3f.gt.0)epsilon=0.7 ! empirical calculation methods for important points ! SK
    if(char3r+snd.le.0 .and. char3f+sndf.gt.0)epsilon=0.7
    if(char3r-snd.le.0 .and. char3f-sndf.gt.0)epsilon=0.7
! LEFT CELL
	epsilon1=epsilon2
	if(char3b.le.0 .and. char3l.gt.0)epsilon1=0.7
    if(char3b+sndb.le.0 .and. char3l+snd.gt.0)epsilon1=0.7
    if(char3b-sndb.le.0 .and. char3l-snd.gt.0)epsilon1=0.7
	
! END CENTRED FAN -- THIS SHOULD BE PLACED BEFORE ANY SUBSEQUENT REDEFINITION OF CHAR3F, ETC
	epsilon=epsilon2
	epsilon1=epsilon2

	q=char3r-gf*p**gamy
	qf=char3f-gf*pf**gamy
	qcf=char3cf-gf*pcf**gamy
	qcsf=char3csf-gf*pcsf**gamy

	char2f=char3csf-soundf	
	qwellf=(qcsf-qcf)/(0.5d0*dt)+char2f*(qf-q)/ds2
	qmaxf=dmax1(q,qcf,qf)+dt*qwellf
	qminf=dmin1(q,qcf,qf)+dt*qwellf

	qnl=(2d0*qcsf-qf*(1.0d0-epsilon))/(1.0d0+epsilon)
	if(qnl>qmaxf) qnl=qmaxf
	if(qnl<qminf) qnl=qminf
!SK  + 4 TERMS IN THE OTHER DIRECTION  (B)
        q=char3r-gb*p**gamy      
        qb=char3b-gb*pb**gamy
        qcb=char3cb-gb*pcb**gamy
        qcsb=char3csb-gb*pcsb**gamy
! SK + EXTRAPOLATION AND CORRECTION IN THE OTHER DIRECTION (B)

    char2b=char3csb-soundb	
	qwellb=(qcsb-qcb)/(0.5d0*dt)+char2b*(q-qb)/ds1
	qmaxb=dmax1(q,qcb,qb)+dt*qwellb
	qminb=dmin1(q,qcb,qb)+dt*qwellb
	qnr=(2d0*qcsb-qb*(1.0d0-epsilon1))/(1.0d0+epsilon1)
	if(qnr>qmaxb) qnr=qmaxb
	if(qnr<qminb) qnr=qminb

! SK

	r=char3l+gb*p**gamy
	rb=char3b+gb*pb**gamy
	rcb=char3cb+gb*pcb**gamy
	rcsb=char3csb+gb*pcsb**gamy

	char1b=char3csb+soundb
	rwellb=(rcsb-rcb)/(0.5d0*dt)+char1b*(r-rb)/ds1	
	rmaxb=dmax1(r,rcb,rb)+dt*rwellb
	rminb=dmin1(r,rcb,rb)+dt*rwellb

        rnr=(2d0*rcsb-rb*(1.0d0-epsilon1))/(1.0d0+epsilon1)
	if(rnr>rmaxb) rnr=rmaxb
	if(rnr<rminb) rnr=rminb
!SK +4 TERMS IN THE OTHER DIRECTION (F)
        r=char3r+gf*p**gamy
        rf=char3f+gf*pf**gamy
        rcf=char3cf+gf*pcf**gamy
        rcsf=char3csf+gf*pcsf**gamy
! SK -- THE SAME FOR THE OTHER DIRECTION (F)

        char1f=char3csf+soundf
	rwellf=(rcsf-rcf)/(0.5d0*dt)+char1f*(rf-r)/ds2	
	rmaxf=dmax1(r,rcf,rf)+dt*rwellf
	rminf=dmin1(r,rcf,rf)+dt*rwellf
        rnl=(2d0*rcsf-rf*(1.0d0-epsilon))/(1.0d0+epsilon)
	if(rnl>rmaxf) rnl=rmaxf
	if(rnl<rminf) rnl=rminf


! SK 
!    if (k.eq.5051) write(14,*) 'u=',u
!    if (k.eq.5051) write(14,*) 'ub=',ub
!    if (k.eq.5051) write(14,*) 'uf=',uf
!    if (k.eq.5051) write(14,*) 'ucb=',ucb
!    if (k.eq.5051) write(14,*) 'ucf=',ucf
!    if (k.eq.5051) write(14,*) 'ucsb=',ucsb
!    if (k.eq.5051) write(14,*) 'ucsf=',ucsf
!	if (k.eq.5051) write(14,*) 'rwell=',r,rb,q,qf,char2b,char2f



    ! SK ALL CHAR'S AT THE RHS BELOW MUST BE REFERRED TO THE CENTRES 
    sound=soundf+soundb
    charc3=char3csf+char3csb
    charc1=(char3l+char3r)/2d0+sound   
    charc2=(char3l+char3r)/2d0-sound



    ! SK DEFINITION //////////////////////////////////////
    gr=gb
    gl=gf
    char3f=char3csf
    char3b=char3csb
    ! END DEFINITION//////////////////////////////////////

    ! SK THE BASE CASE: SUBSONIC LEFT/RIGHT GOING WAVES//////////////////////

    gf=gl
    gb=gr
    rn=rnr
    qn=qnl


! SK MINIMUM VARIATION -- CHAR1


    if(char1f<0.0 .and. char1b>=0.0) then
        rc1=rcb-dt*char1b*rnl/ds1
        rc2=rcf+dt*char1f*rnl/ds2
        var1=dabs(rc1-rc2)

        rc3=rcb-dt*char2b*rnr/ds1
        rc4=rcf+dt*char2f*rnr/ds2
        var2=dabs(rc3-rc4)

        if(var1.lt.var2) then
            rn=rnl
            gb=gl
        else
            rn=rnr
            gb=gr
        end if
    end if

 ! SK END OF MINIMUM VAR CHAR1 

! SK MINIMUM VARIATION -- CHAR2


    if(char2f<0.0 .and. char2b>=0.0) then
        !	QCSLN1=QCL-DT*QLN*CHAR2L/SMESHL*DLNNL
        !	QCSN1=QC+DT*QLN*CHAR2R/SMESH*DLNN
        qc1=qcb-dt*char2b*qnl/ds1
        qc2=qcf+dt*char2f*qnl/ds2
        var1=dabs(qc1-qc2)

        !   QCSLN2=QCL-DT*QRN*CHAR2L/SMESHL*DLNNL
        !	QCSN2=QC+DT*QRN*CHAR2R/SMESH*DLNN
        qc3=qcb-dt*char2b*qnr/ds1
        qc4=qcf+dt*char2f*qnr/ds2
        var2=dabs(qc3-qc4)

        if(var1.lt.var2) then
            qn=qnl
            gb=gl
        else
            qn=qnr
            gb=gr
        end if
    end if

 ! SK END OF MINIMUM VAR CHAR2 

    u1x=qn+gf*(rn-qn)/(gf+gb)
    call GetAlfa(u1x,alfa)
    if(alfa>=1.0d0) then   ! if(charc3<=0.0) then ! SK changed, 2011-05-17
        sn=deff(s,sf,scf,scsf,char3csf,dt,ds2,epsilon)
        uq(1,9)=0.0d0
         uq(1,9)=deff(uq(1,5),uq(1,6),uq(1,7),uq(1,8),char3csf,dt,ds2,epsilon)
       uq(2,9)=deff(uq(2,5),uq(2,6),uq(2,7),uq(2,8),char3csf,dt,ds2,epsilon)
        uq(3,9)=deff(uq(3,5),uq(3,6),uq(3,7),uq(3,8),char3csf,dt,ds2,epsilon)
        pn=((rn-qn)/(gb+gf)/l1r)**(2.0d0*gam/(gam-1.0))
        ron=(pn/exp(sn))**(1.0d0/gam)
	call invg(ar,gi)
		iback=0

    else if(alfa<=0.0d0) then
        sn=defb(s,sb,scb,scsb,char3csb,dt,ds1,epsilon1)
        uq(1,9)=0.0d0
        uq(1,9)=defb(uq(1,1),uq(1,2),uq(1,3),uq(1,4),char3csb,dt,ds1,epsilon1)
        uq(2,9)=defb(uq(2,1),uq(2,2),uq(2,3),uq(2,4),char3csb,dt,ds1,epsilon1)
        uq(3,9)=defb(uq(3,1),uq(3,2),uq(3,3),uq(3,4),char3csb,dt,ds1,epsilon1)
        pn=((rn-qn)/(gb+gf)/l1l)**(2.0d0*gam/(gam-1.0d0))
        ron=(pn/exp(sn))**(1d0/gam)
 	call invg(al,gi)
	else
        sn1 = deff(s,sf,scf,scsf,char3csf,dt,ds2,epsilon)
	ut11 = deff(uq(1,5),uq(1,6),uq(1,7),uq(1,8),char3csf,dt,ds2,epsilon)
        ut12 = deff(uq(2,5),uq(2,6),uq(2,7),uq(2,8),char3csf,dt,ds2,epsilon)
        ut13 = deff(uq(3,5),uq(3,6),uq(3,7),uq(3,8),char3csf,dt,ds2,epsilon)
        pn1 = ((rn-qn)/(gb+gf))**(2.0d0*gam/(gam-1.0))
        ron1 = (pn1/exp(sn1))**(1.0d0/gam)

        sn2 = defb(s,sb,scb,scsb,char3csb,dt,ds1,epsilon1)
	ut21 = defb(uq(1,1),uq(1,2),uq(1,3),uq(1,4),char3csb,dt,ds1,epsilon1)
        ut22 = defb(uq(2,1),uq(2,2),uq(2,3),uq(2,4),char3csb,dt,ds1,epsilon1)
        ut23 = defb(uq(3,1),uq(3,2),uq(3,3),uq(3,4),char3csb,dt,ds1,epsilon1)
        pn2 = ((rn-qn)/(gb+gf))**(2.0d0*gam/(gam-1.0d0))
        ron2 = (pn2/exp(sn2))**(1d0/gam)

        sn = alfa * sn1 + (1.0d0 - alfa) * sn2
        uq(1,9)=0.0d0
        uq(1,9) = alfa * ut11 + (1.0d0 - alfa) * ut21
        uq(2,9) = alfa * ut12 + (1.0d0 - alfa) * ut22
        uq(3,9) = alfa * ut13 + (1.0d0 - alfa) * ut23
        pn = alfa * pn1 + (1.0d0 - alfa) * pn2
        ron = alfa * ron1 + (1.0d0 - alfa) * ron2
		!!!!! NO next
		call invg(ar,gi)
		iback=0
    end if

! SK END OF SUBSONIC ////////////////////////////////////////////

! SK RIGHT-GOING SUPERSONIC/////////////////////////////////////////////


    if(char1b>0.0 .and. char2f>0.0) then
        gb=gr
        gf=gb
        rn=rnr
        qn=qnr

        u1x=0.5d0*(rn+qn)
        sn=defb(s,sb,scb,scsb,char3csb,dt,ds1,epsilon1)
        uq(1,9)=0.0
        uq(1,9)=defb(uq(1,1),uq(1,2),uq(1,3),uq(1,4),char3csb,dt,ds1,epsilon1)
        uq(2,9)=defb(uq(2,1),uq(2,2),uq(2,3),uq(2,4),char3csb,dt,ds1,epsilon1)
        uq(3,9)=defb(uq(2,1),uq(2,2),uq(2,3),uq(2,4),char3csb,dt,ds1,epsilon1)
        pn=((rn-qn)/2.0d0/gb/l1l)**(2.0d0*gam/(gam-1.0d0))
        ron=(pn/exp(sn))**(1d0/gam)
	call invg(al,gi)
	end if


! SK END OF R/G SUPERSONIC ////////////////////////////////

! SK LEFT-GOING SUPERSONIC/////////////////////////////////////////////

    if(char1b<=0.0 .and. char2f<0.0) then
        gb=gl
        gf=gb
        rn=rnl
        qn=qnl

        u1x=0.5d0*(rn+qn)
                sn=deff(s,sf,scf,scsf,char3csf,dt,ds2,epsilon)
	        uq(1,9)=0.0
        uq(1,9)=deff(uq(1,5),uq(1,6),uq(1,7),uq(1,8),char3csf,dt,ds2,epsilon)
	        uq(2,9)=deff(uq(2,5),uq(2,6),uq(2,7),uq(2,8),char3csf,dt,ds2,epsilon)
	        uq(3,9)=deff(uq(3,5),uq(3,6),uq(3,7),uq(3,8),char3csf,dt,ds2,epsilon)
        pn=((rn-qn)/2.0d0/gf/l1r)**(2.0d0*gam/(gam-1.0d0))
        ron=(pn/exp(sn))**(1.0d0/gam)
	call invg(ar,gi)
	iback=0
		end if

! SK END OF L/G SUPERSONIC ////////////////////////////////

! SK SPECIAL TREATMENT FOR THE SONIC POINT -- RIGHT-GOING WAVE

    if(char1b>0.0 .and. char2b<=0.0 .and. char2f>=0.0) then

!SMACHL=WICL/SOUNDL
!SMACH=WIC/SOUND
!SM=0.5*(SMACHL+SMACH)
!G=DSQRT(GAM*DEXP(SRN)**(1./GAM))
!PNN=(RRN/(AL+SM*G))**(2.*GAM/(GAM-1.))
!UNN=SM*G*RRN/(AL+SM*G)
!RONN=(PNN/DEXP(SRN))**(1./GAM)

        gb=gr
        gf=gb
        rn=rnr
        qn=qnr
        sn=defb(s,sb,scb,scsb,char3csb,dt,ds1,epsilon1)
        uq(1,9)=0.0
        uq(1,9)=defb(uq(1,1),uq(1,2),uq(1,3),uq(1,4),char3csb,dt,ds1,epsilon1)
        uq(2,9)=defb(uq(2,1),uq(2,2),uq(2,3),uq(2,4),char3csb,dt,ds1,epsilon1)
        uq(3,9)=defb(uq(3,1),uq(3,2),uq(3,3),uq(3,4),char3csb,dt,ds1,epsilon1)
        gn=dsqrt(gam*dexp(sn)**(1./gam))

        smachl=char3b/soundb
        smachr=char3f/soundf
        smach=0.5*(smachl+smachr)

        pn=(rn/(gb/l1l+gn/l1r*smach))**(2.0d0*gam/(gam-1.0))
        u1x=smach*gn*rn/(gb+gn*smach)
        ron=(pn/exp(sn))**(1.0d0/gam)

	call invg(al,gi)
    end if

! SK END OF SONIC POINT -- R/G WAVE ////////////////////////////////

! SK SPECIAL TREATMENT FOR THE SONIC POINT -- LEFT-GOING WAVE

    if(char1b<=0.0 .and. char1f>=0.0 .and. char2f<0.0) then

        gb=gl
        gf=gb
        rn=rnl
        qn=qnl
        sn=deff(s,sf,scf,scsf,char3csf,dt,ds2,epsilon)
        uq(1,9)=0.0d0
        uq(1,9)=deff(uq(1,5),uq(1,6),uq(1,7),uq(1,8),char3csf,dt,ds2,epsilon)
	        uq(2,9)=deff(uq(2,5),uq(2,6),uq(2,7),uq(2,8),char3csf,dt,ds2,epsilon)
	        uq(3,9)=deff(uq(3,5),uq(3,6),uq(3,7),uq(3,8),char3csf,dt,ds2,epsilon)
        gn=dsqrt(gam*dexp(sn)**(1./gam))

        smachl=char3b/soundb
        smachr=char3f/soundf
        smach=0.5*(smachl+smachr)

        !G=DSQRT(GAM*DEXP(SLN)**(1./GAM))
        !PNN=(-QLN/(AR-SM*G))**(2.*GAM/(GAM-1.))
        !UNN=-SM*G*QLN/(AR-SM*G)

        pn=(-qn/(gb/l1l-gn/l1r*smach))**(2.0d0*gam/(gam-1.0))
        u1x=-smach*gn*qn/(gb-gn*smach)
        ron=(pn/exp(sn))**(1.0d0/gam)
	call invg(ar,gi)
	iback=0
    end if

! SK END OF SONIC POINT -- L/G WAVE ////////////////////////////////

	un=0d0
	if (iback.eq.1) then
	uq(2,9)=uq(2,9)+teta2Alfa1l/teta1Alfa1l*uq(1,9)
	uq(3,9)=uq(3,9)+teta3Alfa1l/teta1Alfa1l*uq(1,9)
	uq(1,9)=u1x
        do i=1,3
	do j=1,3    
     un(i)=un(i)+gi(j,i)*uq(j,9)
	end do
	end do
	else
	uq(2,9)=uq(2,9)+teta2Alfa1r/teta1Alfa1r*uq(1,9)
	uq(3,9)=uq(3,9)+teta3Alfa1r/teta1Alfa1r*uq(1,9)
	uq(1,9)=u1x
        do i=1,3
	do j=1,3    
     un(i)=un(i)+gi(j,i)*uq(j,9)
	end do
	end do
	end if

	vtmp=dsqrt(un(1)**2+un(2)**2+un(3)**2)
    sound=dsqrt(gam*pn/ron)
    dtmint=min(ds1,ds2)/(sound+vtmp)
	if (dtmint.lt.dtmin) dtmin=dtmint
    	if (k.eq.5051) then
!write(14,'(a4,i5,100d12.5)') 'rank',rank,es,char3csb+char3csf
!write(14,'(a4,i5,100d12.5)') 'rank',rank,qn,rn,u1x
    	end if




    return 
end


subroutine usdefnc( un, u,  ub, uf, ucb, ucsb, ucf, ucsf &
                         , pn,p, pb, pf, pcb, pcsb, pcf, pcsf   &
                         , ron,rho,rhob,rhof,rhocb,rhocsb,rhocf,rhocsf   &
                         ,ds1,ds2,gl,gr &
!e1l,e2l,e3l,e1r,e2r,e3r
,dt,dtmin,epsilon2)
implicit none
	!IMPLICIT REAL *8(a-h,o-z)
	integer k,rank
!    real *8 e1l(3),e2l(3),e3l(3),e1r(3),e2r(3),e3r(3)
	real *8 gl(4,3),gr(4,3)
! gl(1,1:3)=e1,sx=gl(1,1)*gl(1,4)
	real *8  un(3), u(3),  ub(3), uf(3), ucb(3), ucsb(3), ucf(3), ucsf(3)
	real *8 pn,p, pb, pf, pcb, pcsb, pcf, pcsf
	real *8 rhon,rho,rhob,rhof,rhocb,rhocsb,rhocf,rhocsf
	real *8 ds1,ds2,dt,dtmin
	real *8 ut(3,10)
    real *8 ei,z1,vtmp,dttmp
    real *8 char3,q,r,qn,rn,u1x,ron,vf,wf,vb,wb
    real *8 char3f,soundf,gf,sf,qf,qcf,qcsf,char2f,qwellf,qmaxf,qminf,qwell
    real *8 char3cf
    real *8 char3csf
    real *8 char3b,soundb,gb,sb,rb,rcb,rcsb,char2b,rmaxb,rminb
    real *8 char3cb
    real *8 char3csb
! SK EXTRA VARIABLES
    real *8 qb,qcb,qcsb,rf,rcf,rcsf,rwellb
	real *8 qwellb,qmaxb,qminb,rwellf,rmaxf,rminf,rwell
	real *8 rnr,rnl,qnl,qnr,char1f,char1b,snd,sndf,sndb
	real *8 smachl,smachr,smach,gn,rc1,rc2,rc3,rc4,qc1,qc2,qc3,qc4,var1,var2
! END SK
    real *8 charc1,charc2,charc3
    real *8 a(3,3)
    real *8 epsilon
    real *8 es(3),gi(3,3)
    real *8 deff,defb
    real *8 dtmint,gam,gam2,gamx,gamy,s,scb,scf,scsb,scsf,sn,sound
    integer i,j
    external deff
    external defb
    real *8 epsilon1,epsilon2
real *8	teta1Alfa1l,l1l
real *8	teta2Alfa1l
real *8	teta3Alfa1l
real *8	teta1Alfa1r,l1r
real *8	teta2Alfa1r
real *8	teta3Alfa1r


    real(8) :: alfa, sn1,sn2, ut12,ut13,ut22,ut23, pn1,pn2, ron1,ron2,char3l,char3r
    epsilon=epsilon2
    gam=1.4d0
    gamx=2.0d0*dsqrt(gam)/(gam-1.0d0)
    gamy=((gam-1.0d0)/(2.0d0*gam))
    gam2=2.0d0*gam
                               ! first index - ort
	call invg(gl,gi)
	teta1Alfa1l=(gi(1,1)*gl(1,1)+gi(2,1)*gl(1,2)+gi(3,1)*gl(1,3))*gl(4,1)
	teta2Alfa1l=(gi(1,1)*gl(2,1)+gi(2,1)*gl(2,2)+gi(3,1)*gl(2,3))*gl(4,1)
	teta3Alfa1l=(gi(1,1)*gl(3,1)+gi(2,1)*gl(3,2)+gi(3,1)*gl(3,3))*gl(4,1)
	l1l=dsqrt(teta1Alfa1l/gl(4,1))
	ut=0d0
! 1-u 2-ub  4-ucb 5-ucsb
	do i=1,3
	do j=1,3
	ut(i,1)=ut(i,1)+gl(j,i)*u(j)
	ut(i,2)=ut(i,2)+gl(j,i)*ub(j)
	ut(i,3)=ut(i,3)+gl(j,i)*ucb(j)
	ut(i,4)=ut(i,4)+gl(j,i)*ucsb(j)
	end do
	end do
	do j=1,4
	ut(2,j)=ut(2,j)-teta2Alfa1l/teta1Alfa1l*ut(1,1)
	ut(3,j)=ut(3,j)-teta3Alfa1l/teta1Alfa1l*ut(1,1)
	end do
    char3l =ut(1,1)
    char3b=ut(1,2)
    char3cb=ut(1,3) 
    char3csb=ut(1,4)
s=dlog(p)-gam*dlog(rho)
sb=dlog(pb)-gam*dlog(rhob)
scb=dlog(pcb)-gam*dlog(rhocb)
scsb=dlog(pcsb)-gam*dlog(rhocsb)
	gb=gamx*dexp(scsb/gam2)*l1l
soundb=dsqrt(gam*pcsb/rhocsb)
	r=char3l+gb*p**gamy
	rb=char3b+gb*pb**gamy
	rcb=char3cb+gb*pcb**gamy
	rcsb=char3csb+gb*pcsb**gamy
	char2b=char3csb+soundb
	rwell=(rcsb-rcb)/(0.5d0*dt)+char2b*(r-rb)/ds1
	
	rmaxb=dmax1(r,rcb,rb)+dt*rwell
	rminb=dmin1(r,rcb,rb)+dt*rwell
        rn=(2d0*rcsb-rb*(1.0d0-epsilon))/(1.0d0+epsilon)
	if(rn>rmaxb) rn=rmaxb
	if(rn<rminb) rn=rminb
	call invg(gr,gi)
!	teta1Alfa1r= (gr(1,1)*gr(1,1)+gr(2,1)*gr(1,2)+gr(3,1)*gr(1,3))*gr(4,1)
!	teta2Alfa1r= (gr(1,1)*gr(2,1)+gr(2,1)*gr(2,2)+gr(3,1)*gr(2,3))*gr(4,1)
!	teta3Alfa1r= (gr(1,1)*gr(3,1)+gr(2,1)*gr(3,2)+gr(3,1)*gr(3,3))*gr(4,1)

	teta1Alfa1r= (gi(1,1)*gr(1,1)+gi(2,1)*gr(1,2)+gi(3,1)*gr(1,3))*gr(4,1)
	teta2Alfa1r= (gi(1,1)*gr(2,1)+gi(2,1)*gr(2,2)+gi(3,1)*gr(2,3))*gr(4,1)
	teta3Alfa1r= (gi(1,1)*gr(3,1)+gi(2,1)*gr(3,2)+gi(3,1)*gr(3,3))*gr(4,1)


	l1r=dsqrt(teta1Alfa1r/gr(4,1))
	
! 6-u 7-uf  9-ucf 10-ucsf
	do i=1,3
	do j=1,3
	ut(i,5)=ut(i,5)+gr(j,i)*u(j)
	ut(i,6)=ut(i,6)+gr(j,i)*uf(j)
	ut(i,7)=ut(i,7)+gr(j,i)*ucf(j)
	ut(i,8)=ut(i,8)+gr(j,i)*ucsf(j)
	end do
	end do
	do j=5,8
	ut(2,j)=ut(2,j)-teta2Alfa1r/teta1Alfa1r*ut(1,5)
	ut(3,j)=ut(3,j)-teta3Alfa1r/teta1Alfa1r*ut(1,5)
	end do

    char3r =ut(1,5)
    char3f=ut(1,6)
    char3cf=ut(1,7) 
    char3csf=ut(1,8)
s=dlog(p)-gam*dlog(rho)
sf=dlog(pf)-gam*dlog(rhof)
scf=dlog(pcf)-gam*dlog(rhocf)
scsf=dlog(pcsf)-gam*dlog(rhocsf)
	gf=gamx*dexp(scsf/gam2)*l1r
soundf=dsqrt(gam*pcsf/rhocsf)

	q=char3r-gf*p**gamy
	qf=char3f-gf*pf**gamy
	qcf=char3cf-gf*pcf**gamy
	qcsf=char3csf-gf*pcsf**gamy
	char2f=char3csf-soundf	
	qwell=(qcsf-qcf)/(0.5d0*dt)+char2f*(qf-q)/ds2
	qmaxf=dmax1(q,qcf,qf)+dt*qwell
	qminf=dmin1(q,qcf,qf)+dt*qwell
!	qn=2*qcsf-qf
	qn=(2d0*qcsf-qf*(1.0d0-epsilon))/(1.0d0+epsilon)
	if(qn>qmaxf) qn=qmaxf
	if(qn<qminf) qn=qminf


sound=soundf+soundb
charc3=char3f+char3b
charc1=charc3+sound   
charc2=charc3-sound
if(charc1>0.0 .and. charc2>0.0) then
u1x=0.5d0*(rn+qn)
 sn=defb(s,sb,scb,scsb,char3csb,dt,ds1,epsilon)
	ut(1,9)=0.0
	ut(2,9)=defb(ut(2,1),ut(2,2),ut(2,3),ut(2,4),char3csb,dt,ds1,epsilon)
	ut(3,9)=defb(ut(3,1),ut(3,2),ut(3,3),ut(3,4),char3csb,dt,ds1,epsilon)
	ut(2,9)=ut(2,9)+teta2Alfa1l/teta1Alfa1l*u1x
	ut(3,9)=ut(3,9)+teta3Alfa1l/teta1Alfa1l*u1x
	call invg(gl,gi)
pn=((rn-qn)/2.0d0/gb/l1l)**(2.0d0*gam/(gam-1.0d0))
ron=(pn/exp(sn))**(1d0/gam)
end if
if(charc1<=0.0 .and. charc2<0.0) then
u1x=0.5d0*(rn+qn)
        sn=deff(s,sf,scf,scsf,char3csf,dt,ds2,epsilon)
	ut(1,9)=0.0
	ut(2,9)=deff(ut(2,5),ut(2,6),ut(2,7),ut(2,8),char3csf,dt,ds2,epsilon)
	ut(3,9)=deff(ut(3,5),ut(3,6),ut(3,7),ut(3,8),char3csf,dt,ds2,epsilon)
	ut(2,9)=ut(2,9)+teta2Alfa1r/teta1Alfa1r*u1x
	ut(3,9)=ut(3,9)+teta3Alfa1r/teta1Alfa1r*u1x
	call invg(gr,gi)
pn=((rn-qn)/2.0d0/gf/l1r)**(2.0d0*gam/(gam-1.0d0))
ron=(pn/exp(sn))**(1.0d0/gam)
end if
!if(charc1>0.0 .and. charc2<=0.0) then
u1x=qn+gf*(rn-qn)/(gf+gb)
if(char3csb+char3csf<=0.0) then
u1x=qn+gf*(rn-qn)/(gf+gb)
        sn=deff(s,sf,scf,scsf,char3csf,dt,ds2,epsilon)
	ut(1,9)=0.0
	ut(2,9)=deff(ut(2,5),ut(2,6),ut(2,7),ut(2,8),char3csf,dt,ds2,epsilon)
	ut(3,9)=deff(ut(3,5),ut(3,6),ut(3,7),ut(3,8),char3csf,dt,ds2,epsilon)
	ut(2,9)=ut(2,9)+teta2Alfa1r/teta1Alfa1r*u1x
	ut(3,9)=ut(3,9)+teta3Alfa1r/teta1Alfa1r*u1x
	call invg(gr,gi)
pn=((rn-qn)/(gb+gf)/l1r)**(2.0d0*gam/(gam-1.0))
ron=(pn/exp(sn))**(1.0d0/gam)
else 
u1x=qn+gf*(rn-qn)/(gf+gb)
 sn=defb(s,sb,scb,scsb,char3csb,dt,ds1,epsilon)
	ut(1,9)=0.0
	ut(2,9)=defb(ut(2,1),ut(2,2),ut(2,3),ut(2,4),char3csb,dt,ds1,epsilon)
	ut(3,9)=defb(ut(3,1),ut(3,2),ut(3,3),ut(3,4),char3csb,dt,ds1,epsilon)
	ut(2,9)=ut(2,9)+teta2Alfa1l/teta1Alfa1l*u1x
	ut(3,9)=ut(3,9)+teta3Alfa1l/teta1Alfa1l*u1x
	call invg(gl,gi)
pn=((rn-qn)/(gb+gf)/l1r)**(2.0d0*gam/(gam-1.0d0))
ron=(pn/exp(sn))**(1d0/gam)
end if
	un=0d0
	ut(1,9)=u1x

        do i=1,3
	do j=1,3    
     un(i)=un(i)+gi(i,j)*ut(j,9)
	end do
	end do



	return
	end
	subroutine invg(g,a)
	real *8 g(4,3),a(3,3)
!       a(1,1) a(1,2),a(1,3)
!       a(2,1) a(2,2),a(2,3)
!       a(3,1) a(3,2),a(3,3)
      REAL *8 B(9),C1
	do j=1,3
	b(j)=g(1,j)
	b(j+3)=g(2,j)
	b(j+6)=g(3,j)
	end do
!
      C1=(B(5)*B(9)-B(6)*B(8))*B(1)+(B(6)*B(7)-B(4)*B(9))*B(2)+(B(4)*B(8)-B(5)*B(7))*B(3)
!      A(1,1)=(B(5)*B(9)-B(6)*B(8))/C1
!      A(1,2)=(B(6)*B(7)-B(4)*B(9))/C1
!      A(1,3)=(B(4)*B(8)-B(5)*B(7))/C1
!      A(2,1)=(B(8)*B(3)-B(9)*B(2))/C1
!      A(2,2)=(B(9)*B(1)-B(7)*B(3))/C1
!      A(2,3)=(B(7)*B(2)-B(8)*B(1))/C1
!      A(3,1)=(B(2)*B(6)-B(3)*B(5))/C1
!      A(3,2)=(B(3)*B(4)-B(1)*B(6))/C1
!      A(3,3)=(B(1)*B(5)-B(2)*B(4))/C1
      A(1,1)=(B(5)*B(9)-B(6)*B(8))/C1
      A(2,1)=(B(6)*B(7)-B(4)*B(9))/C1
      A(3,1)=(B(4)*B(8)-B(5)*B(7))/C1
      A(1,2)=(B(8)*B(3)-B(9)*B(2))/C1
      A(2,2)=(B(9)*B(1)-B(7)*B(3))/C1
      A(3,2)=(B(7)*B(2)-B(8)*B(1))/C1
      A(1,3)=(B(2)*B(6)-B(3)*B(5))/C1
      A(2,3)=(B(3)*B(4)-B(1)*B(6))/C1
      A(3,3)=(B(1)*B(5)-B(2)*B(4))/C1

!
	return
	end

	subroutine smothg(x,y)
	real *8 x(4,3),y(4,3)
	real *8 c1,c2,ct,r(4),z(4,3)
	do j=1,3
	x(4,j)=dsqrt(x(1,j)**2+x(2,j)**2+x(3,j)**2)
	x(1:3,j)=x(1:3,j)/x(4,j)
	end do
	do j=1,3
	y(4,j)=dsqrt(y(1,j)**2+y(2,j)**2+y(3,j)**2)
	y(1:3,j)=y(1:3,j)/y(4,j)
	end do
	c1=x(1,2)*y(1,2)+x(2,2)*y(2,2)+x(3,2)*y(3,2)
	c2=x(1,2)*y(1,3)+x(2,2)*y(2,3)+x(3,2)*y(3,3)
	if (dabs(c1).lt.dabs(c2)) then
	r(1:4)=y(1:4,2)
	y(1:4,2)=y(1:4,3)
	y(1:4,3)=r(1:4)
	ct=c1
	c1=c2
	c2=ct
	end if
	if (c1.lt.0d0) y(1:3,2)=-y(1:3,2)
	c1=x(1,3)*y(1,3)+x(2,3)*y(2,3)+x(3,3)*y(3,3)
	if (c1.lt.0d0) y(1:3,3)=-y(1:3,3)
	z=(x+y)/2d0
	do j=1,3
	c1=dsqrt(z(1,j)**2+z(2,j)**2+z(3,j)**2)
	z(1:3,j)=z(1:3,j)/c1

	end do
	x=z
	y=z
	return
	end
